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Progress  Report  on  the  MR  Elastography  System 
for  Breast  Cancer  Detection 


1.0  INTRODUCTION 

Breast  cancer  is  the  most  commonly  diagnosed  cancer  in  women.  Early  diagnosis,  which 
is  critical  for  favorable  clinical  outcomes,  is  complicated  by  the  difficulty  of  detecting  small 
tumors.  One  physical  property  that  clearly  distinguishes  healthy  from  cancerous  tissue  is 
mechanical  stiffness  (hardness).  For  this  reason,  palpation  has  long  been  used  for  early 
screening.  Researchers  have  attempted  to  combine  external  mechanical  stimulation  and 
Magnetic  Resonance  Imaging  (MRI)  to  quantitatively  measure  the  Young’s  modulus  of  tissue 
throughout  both  the  breast  and  the  prostate.  This  technique.  Magnetic  Resonance  Elastography 
(MRE)  has  been  called  “palpation  at  a  distance.”  One  of  the  most  challenging  technical  issues  of 
MRE  is  the  solution  of  the  “inverse  problem,”  i.e.,  quantitatively  inferring  Young’s  modulus 
from  MRI-measured  tissue  displacement  data.  Addressing  this  problem  was  the  focus  of 
Creare’s  efforts  in  the  first  year  of  this  project. 

2.0  BACKGROUND 

Breast  cancer  is  the  most  commonly  diagnosed  cancer  in  women,  with  a  current  mortality 
rate  of  approximately  23  per  100,000.  Early  diagnosis,  which  is  critical  for  favorable  clinical 
outcomes,  is  complicated  by  the  difficulties  associated  with  detecting  small  tumors.  While 
mammography  is  very  useful,  the  X-ray  attenuation  properties  of  tumors  and  normal  tissues  are 
similar,  often  making  it  difficult  to  distinguish  small  tumors  from  healthy  tissue.  The  same 
problems  afflict  ultrasound,  since  the  acoustic  impedance  of  tumors  is  only  slightly  different 
from  that  of  normal  tissue. 

One  physical  characteristic  that  does  clearly  distinguish  healthy  from  cancerous  tissue  is 
mechanical  stiffness  (hardness).  The  Young’s  modulus  of  breast  tumors  can  be  one  to  two 
orders  of  magnitude  larger  than  that  of  normal  tissue  (Muthupillai  and  Ehman,  1996).  For  this 
reason,  palpation  has  long  been  used  by  physicians  and  patients  for  early  screening.  However, 
manual  palpation  is  not  very  effective  for  tumors  lying  deep  within  the  breast  and  is  not 
quantitative.  For  these  reasons,  over  the  past  decade,  numerous  efforts  have  been  underway  to 
combine  external  mechanical  stimulation  and  various  imaging  techniques  to  quantitatively 
measure  the  Young’s  modulus  of  tissue  throughout  both  the  breast  and  the  prostate.  This  has 
been  termed  “palpation  at  a  distance.”  Because  of  its  inherent  sensitivity,  much  current  research 
is  focused  on  using  magnetic  resonance  imaging  (MRI),  and  the  overall  process  is  termed 
magnetic  resonance  elastography  (MRE). 

Probably  the  most  challenging  technical  aspect  of  MRE  is  the  solution  of  the  “inverse 
problem,”  i.e.,  quantitatively  inferring  Young’s  modulus  from  MRI-measured  tissue 
displacement  data  (Trahey,  2001).  Addressing  this  problem  was  the  focus  of  Creare’s  efforts  in 
the  first  year  of  this  project. 

2.1  Previous  Work 

A  number  of  researchers  have  investigated  use  of  displacement  images  to  infer  the 
Young’s  modulus  of  tissue.  Early  work  focused  on  utilizing  ultrasound  images,  but  much 
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current  research  employs  MRI  which  can  accurately  measure  much  smaller  displacements  (less 
than  one  micron)  and  has  higher  resolution  than  ultrasound. 

MRE  techniques  can  be  broadly  distinguished  by  whether  they  employ  static  or  dynamic 
displacement  of  the  tissue  surface.  Static  techniques  measure  the  change  in  tissue  displacement 
resulting  from  a  quasi-static  push  on  the  outside  of  the  breast,  or  a  sequence  of  pushes.  This 
method  is,  in  principle,  simpler  than  dynamic  methods,  but  certain  mathematical  difficulties  arise 
from  the  use  of  static  displacements.  These  difficulties  appear  to  have  led  to  the  virtual 
abandonment  of  this  technique  (Raghavan  and  Yagle,  1994;  Chenevert  et  al.  1998;  Skovoroda 
et  al.  1995). 

Dynamic  MRE  usually  involves  the  generation  of  oscillatory  motions  at  the  surface  of  the 
breast  to  induce  shear  waves  in  the  tissue.  Because  shear  waves  are  strongly  attenuated, 
researchers  have  also  investigated  the  use  of  ultrasound  to  induce  shear  waves  in  deep  tumors  via 
longitudinal/shear  mode  conversion  (Wu  et  al.  2000).  In  either  case,  the  motion  of  the  interior  of 
the  breast  (or  phantom)  is  then  measured  using  MRI,  and  the  elastic  modulus  is  inferred  from 
these  measurements. 

MRE  requires  the  use  of  relatively  low  excitation  frequencies  since  motion  at  more  than 
about  1  kHz  cannot  be  resolved  by  MRI.  Shear  waves  are  used,  because  even  at  these  low 
frequencies  they  have  a  sufficiently  short  wavelength  to  allow  the  resolution  of  small  features 
such  as  tumors.  Longitudinal  waves  have  a  much  faster  sound  speed,  a  correspondingly  larger 
wavelength,  and  thus  little  resolving  power. 

Two  techniques  have  been  used  to  infer  Young’s  modulus  from  shear-wave-induced 
displacement  data.  Early  research  efforts  employed  direct  measurements  of  the  wavelength  of 
shear  waves  at  various  points  in  a  tissue  phantom  to  determine  shear  wave  speed.  From  shear 
wave  speed,  the  modulus  can  be  readily  determined  (Muthupillai  et  al.  1995;  Manduca  et  al. 
1996).  A  difficulty  of  this  direct  technique  is  that  in  real  tissue,  acoustic  waves  scatter  off  tissue 
inhomogeneities  and  the  scattered  waves  interfere,  generally  leading  to  a  very  complex  pattern  of 
displacements  (van  Houten  et  al.  1999).  Such  a  complex  pattern  does  not  lend  itself  to  direct 
measurements  of  wavelength.  To  avoid  this  problem,  several  researchers  have  turned  to  physics- 
based  modeling  of  tissue  motion.  In  this  method,  the  equations  of  motion  of  the  tissue  are 
defined  in  terms  of  the  known  oscillatory  boundary  conditions  and  a  set  of  assumed  tissue 
properties  (i.e..  Young’s  modulus).  These  equations  are  then  solved,  usually  with  finite  element 
techniques,  and  then  the  calculated  displacement  patterns  are  compared  to  those  measured  with 
MRI.  The  assumed  Young’s  moduli  in  the  model  are  then  systematically  altered  and  the  process 
repeated  until  the  predicted  and  measured  displacement  patterns  match  (to  within  some 
tolerance). 

One  of  the  most  active  groups  pursuing  this  approach  is  at  Dartmouth  College’s  Medical 
School  and  Thayer  School  of  Engineering.  They  have  reported  substantial  success,  but  two  key 
problems  stand  out; 
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1.  The  time  required  to  converge  the  inverse  calculations  can  be  measured  in  days,  even  for 
a  relatively  coarse  mesh  on  a  high-performance  workstation.  While  this  has  led  to  the 
development  of  so-called  subzone  techniques  and  utilization  of  parallel  processing  to 
reduce  exeeution  times,  this  continues  to  represent  a  significant  practical  problem  for 
using  MRE  in  a  clinical  setting  (Van  Houten  et  al.  1999;  Van  Houten,  et  al.  2000a). 

2.  While  good  results  have  been  reported  in  numerical  simulations,  the  results  in  real  tissue 
and  in  phantoms  have  been  decidedly  mixed.  This  has  lead  to  recent  efforts  to  model 
transient  effects  (rather  than  assuming  steady-state),  motion  in  three  dimensions  (rather 
than  the  plane  strain  approximation  used  previously),  and  speculation  that  viscoelastic 
behavior  may  need  to  be  modeled  (Van  Houten  et  al.  2000b). 

2.2  Overview  of  Creare  Efforts 

Creare  efforts  to  date  have  focused  primarily  on  reducing  the  time  required  to  solve  the 
inverse  problem.  Further,  our  most  recent  efforts  have  begun  to  address  how  to  increase  the 
fidelity  of  the  reconstruction  process  (i.e.,  the  underlying  finite  element  model).  In  this  section, 
we  will  introduce  these  areas  of  research.  The  next  two  sections  will  describe  our  work  in  more 
detail. 


The  Dartmouth  group  obtains  updated  Young’s  modulus  estimates  by  minimizing  the 
following  figure  of  merit  of  the  results: 

j  =  (1) 

/=l„n 


In  Equation  (1),  u  and  v  denote  the  x-  and  y-direction  displacements  at  each  of  the  n 
nodes  of  the  model  at  some  specified  time  in  their  oscillatory  cycle.  The  superscript  denotes 
displacements  obtained  from  measured  MRI  data  (which  can  come  from  an  actual  MRI  or  from  a 
finite  element  simulation).  The  lack  of  a  superscript  denotes  code-calculated  displacements 
based  on  the  current  estimate  for  Young’s  modulus  at  each  node.  Note  that  in  this  and  the 
following  discussion,  we  assume  that  a  two-dimensional,  plane  strain  approximation  is  adequate, 
although  the  formulation  can  be  readily  generalized  to  a  full  3-D  analysis.  We  will  later  briefly 
discuss  a  modeling  possibility  that  may  alleviate  the  need  to  consider  a  fully  three-dimensional 
analysis. 

The  Dartmouth  group  minimizes  J  using  a  well-known  technique  known  as  the 
Levenberg-Marquardt  algorithm  (Press  et  al.  1986).  This  requires  the  estimation  of  the  Hessian 
matrix  associated  with  J,  which  in  turn  requires  the  calculation  of  terms  such  as: 


3m, 

w, 


i=l,...n;  k=l,...p 


(2) 


where  Ek  denotes  the  Young’s  modulus  of  the  k*  element  in  the  finite  element  model.  If  there 
are  n  nodes  in  the  model,  the  calculation  of  the  nodal  displacements  requires  a  time  proportional 
to  roughly  n^.  However,  the  caleulation  of  the  terms  (2)  requires  a  time  that  is  p  times  longer 
than  this,  which  for  low  order  elements  (for  which  p  is  roughly  equal  to  n)  is  on  the  order  of  n"*. 
This  is  a  significant  disadvantage  as  the  number  of  elements  and  nodes  grows  unless  the 
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calculation  of  these  quantities  makes  convergence  to  the  proper  Young’s  modulus  distribution 
much  more  rapid. 


Since  the  time  required  to  calculate  the  terms  in  the  approximate  Hessian  matrix  is  very 
long,  we  have  investigated  the  use  of  a  quasi-Newton  method,  specifically  the  Boyden-Fletcher- 
Goldfarb-Shanno  (BFGS)  algorithm  (Press  et  al.  1986).  While  this  method  is  likely  to  require 
more  iterations  to  converge  than  Levenberg-Marquardt,  BFGS  requires  the  calculation  of  only 
the  gradient  of  the  figure-of-merit,  i.e.,  the  terms: 


k=l,...p 


(3) 


To  facilitate  use  of  this  technique,  we  developed  an  algorithm  that  very  efficiently  calculates  the 
terms  of  the  gradient.  This  algorithm  is  described  in  Section  2.3. 


When  we  employed  tissue  mechanical  properties  similar  to  those  used  in  the  literature  in 
numerical  experiments,  we  achieved  good  success  in  estimating  the  values  of  the  Young’s 
modulus  from  simulated  displacement  data.  We  subsequently  determined  that  much  of  the  work 
reported  in  the  literature  appears  to  use  a  value  for  Poisson’s  ratio  that  can  severely  distort  the 
results.  We  are  currently  working  to  remedy  this  problem.  These  efforts  are  described  in 
Section  2.4. 


2.3  Improved  Inverse  Solution  Techniques 
2.3.1  Solution  of  the  Tissue  Equations  of  Motion 

To  develop  an  improved  inverse  solution  technique,  it  is  first  necessary  to  develop  a 
numerical  model  for  the  tissue  displacements  given  an  assumed  set  of  tissue  properties.  While  a 
commercial  product  such  as  ANSYS  can  (and  was)  used  for  preliminary  analyses  (see 
Section  2.4),  ultimately  a  customized  computer  code  is  needed  so  that  the  displacement  solution 
can  be  integrated  with  the  solution  of  the  inverse  problem.  We  developed  a  finite  element  code 
in  Fortran  90  for  this  purpose  using  the  techniques  outlined  in  Smith  and  Griffiths  (1998).  One 
design  goal  was  to  construct  a  code  that  would  be  capable  of  handling  relatively  large  problems 
so  that  this  code  could  serve  as  a  test  bed  to  investigate  techniques  for  accelerating  the  solution 
of  very  large  MRE  inverse  problems. 

The  finite  element  code  calculates  the  steady-state  amplitude  of  a  damped,  linear  elastic 
solid  in  plane  strain,  subjected  to  a  sinusoidal  shear  or  longitudinal  displacement  along  one  or 
more  edges.  The  mesh  is  generated  by  a  separate  program  written  for  this  purpose  and  read  from 
a  file.  Various  types  of  elements  can  be  used,  including  triangles  and  quadrilaterals  of  various 
numbers  of  nodes. 

In  terms  of  the  development  of  a  finite-element  model,  the  equations  of  motion  of  the 
nodes  of  the  model  can  be  written: 


KMu  +  fiKM  —  +  MM^  =  F 
dt  dt^ 


(4) 
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where 

u  =  time-dependent  displacement  vector  of  the  nodes 
KM  =  stiffness  matrix  (relates  displacements  to  stresses) 

MM  =  mass  matrix  (defines  inertia) 

P  =  Rayleigh  damping  coefficient  that  defines  the  damping  ratio  in  terms  of  the  stiffness  matrix 

F  =  time-dependent  applied  external  forces 

It  can  be  shown  that  if  the  desired  damping  ratio  is  y,  then  for  an  angular  excitation  frequency  (O 
one  should  use  (Smith  and  Griffiths,  1998): 


P  =  ^  (5) 

m 

It  should  be  noted  that  it  is  also  possible  to  define  the  damping  in  terms  of  the  stiffness  matrix, 
using  a  slightly  different  form  for  the  coefficient. 

For  the  case  where  we  assume  that  a  steady-state  harmonic  excitation  at  angular  rate  (O  is 
employed,  we  define  the  time-dependent  forces  and  motion  in  terms  of  their  amplitudes  via: 

F{t)  =  F/”  (6) 

u{t)  =  (7) 

Substituting  Equations  (6)  and  (7)  into  Equation  (4),  we  obtain: 

{KM  +  iPmKM  -  m^MM  )u„  =  F„  (8) 

Equation  (8)  is  a  time-independent  matrix  equation  for  the  displacement  amplitudes  in  terms  of 
the  amplitudes  of  the  applied  forces  and  an  assumed  set  of  tissue  properties.  It  is  convenient  to 
write  this  in  the  standard  matrix  equation  form: 


Ax -b  (9) 

where  A  is  the  (complex)  matrix  involving  the  terms  in  parenthesis  in  Equation  (8),  b  represents 
the  (real)  amplitude  of  the  forces  supplied  at  the  surface,  and  x  is  the  solution  vector  which  is 
therefore  also  complex.  Physically,  the  fact  that  x  is  complex  arises  from  differences  in  phase 
between  the  displacements  at  different  points  within  the  solid.  The  variations  in  phase  are 
caused  by  the  effects  of  damping. 

Actually,  one  ordinarily  does  not  know  the  forces  applied  at  the  surface  of  the  breast  but 
rather  the  magnitude  of  the  forced  oscillations.  For  this  reason.  Equation  (9)  is  rewritten  in 
practice  so  that  the  surface  nodes  are  forced  to  have  a  known  amplitude  of  vibration.  One  simple 
way  to  do  this  is  to  replace  the  associated  matrix  elements  in  Equation  (9)  as  follows: 

At  =  A*  (for  >  o^oal  to  externally  excited  degrees-of-freedom)  (10) 
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bj^  -  amplitude  of  excitation  of  k*  degree-of-freedom 

The  terms  in  b  corresponding  to  non-extemally  forced  nodes  are  zero. 

A  difficult  aspect  of  this  problem  is  obtaining  a  solution  of  the  matrix  Equation  (8)  or  (9). 
This  solution  is  somewhat  complicated  in  this  case  because  the  large  matrix  A,  while  sparse,  is 
complex,  usually  symmetric,  and  therefore  not  Hermitian  (i.e.,  A^). 

A  direct  solution  scheme  is  available  for  equations  of  the  form  of  (9)  that  takes  advantage 
of  the  symmetry  of  A  (Harwell,  undated).  This  is  a  specialization  of  the  typical  LU- 
decomposition  scheme  ordinarily  used  to  solve  real  matrix  equations  of  the  form  (8).  An 
advantage  of  such  schemes  is  that  once  the  matrix  A  has  been  decomposed,  e.g.,  factored  into 
lower-  and  upper-triangular  matrices  L  and  U,  solutions  for  different  forcing  vectors  b  can  be 
obtained  with  virtually  no  additional  effort.  This  is  useful  for  calculating  the  gradients,  as 
demonstrated  below. 

However,  for  very  large  models,  so-called  “mesh-free”  techniques  are  attractive, 
especially  when  coupled  to  iterative  matrix  solution  techniques  (Smith  and  Griffiths,  1998).  The 
advantage  of  this  scheme  is  that  the  complete  stiffness  and  mass  matrices  MM  and  KM  need 
never  be  assembled  into  memory,  even  in  a  sparse  format.  This  allows  very  large  problems  to  be 
solved  even  when  memory  is  relatively  limited.  Good  results  were  obtained  in  our  case  solving 
Equation  (9)  using  a  variant  of  the  preconditioned  conjugate  gradient  technique  (PCG),  namely 
the  biconjugate  gradient  squared  iterative  method  (Joubert,  et  al.  undated).  For  example,  the 
solution  of  a  20,000  element  grid  with  20,000  nodes  was  obtained  in  less  than  one  and  one-half 
minutes  on  a  400  MHz  personal  computer. 

One  problem  was  encountered  using  the  iterative  solution  technique.  Strangely, 
convergence  of  these  techniques  can  be  adversely  affected  if  one  begins  the  iteration  process 
with  a  solution  that  is  close  to  the  correct  one.  This  is  potentially  a  serious  disadvantage,  since 
when  performing  the  inverse  problem  one  usually  has  an  excellent  approximation  to  the  correct 
solution  and  use  should  be  made  of  this  to  accelerate  convergence.  That  is,  if  only  a  small 
change  in  Young’s  modulus  is  made,  the  previously  converged  results  will  be  very  close  to  the 
new  ones.  This  problem  was  solved  by  calculating  the  change  in  displacements  from  the 
previous  solution.  Let  A  be  the  matrix  for  which  a  converged  solution  x  has  been  obtained  for  a 
constant  forcing  function  b.  Let  6A  be  the  change  caused  by  altered  assumptions  on  Young’s 
modulus.  The  change  in  the  solution,  5x,  is  given  by: 

{A  +  dA){x  +  dx)  =  b  (12) 

The  previous  converged  solution  satisfied  Equation  (9).  Subtracting  Equation  (9)  from 
Equation  (12)  and  keeping  only  the  first  order  terms,  we  have: 

Adx  =  -5Ax  (13) 

The  advantage  of  this  formulation  is  that  the  iteration  process  can  begin  with  =  0 ,  which  is 
required  for  good  convergence  with  the  PCG  method,  but  the  process  completes  much  quicker 
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since  the  required  accuracy  in  &c  is  much  smaller  than  is  needed  for  the  complete  solution 


An  exact  solution  of  the  displacements  is  available  for  the  special  case  of  a  rectangular 
solid,  fixed  on  all  but  one  edge  and  subjected  to  a  transverse  excitation  whose  amplitude  varies 
sinusoidally  across  the  unfixed  edge  (Gao,  1995).  This  solution  was  compared  to  the  predictions 
of  the  finite  element  code.  The  latter  involved  a  10  cm  x  10  cm  square,  modeled  using  1  mm 
X  1  mm  3-node  triangular  elements.  Figure  1  shows  the  displacement  amplitude  parallel  to  the 
direction  of  excitation  down  the  midline  of  the  solid.  The  code  results  agree  reasonably  well 
with  the  exact  solution.  It  is  noteworthy  that  the  value  of  the  damping  ratio  used  in  this  example 
problem  (^  =  0.1)  is  not  quite  sufficient  to  prevent  some  standing  waves  from  developing  along 
the  side  of  the  rectangle  that  is  opposite  the  driven  edge. 

Of  course,  the  accuracy  of  the  solution  is  controlled  by  the  fineness  of  the  grid  and  the 
fidelity  afforded  by  the  various  types  of  elements.  Figure  1  utilized  a  fairly  coarse  mesh,  and 
more  accurate  solutions  can  be  obtained  by  using  more  elements,  quadrilateral  elements  instead 
of  3-node  triangles,  or  triangles  with  more  than  3  nodes.  Figure  2  shows  the  behavior  of  a 
snapshot  of  the  displacements  at  a  particular  time,  i.e.,  the  real  part  of  the  complex  displacements 
(not  their  total  amplitude  as  shown  in  Figure  1).  The  cases  with  a  finer  mesh  (i.e.,  more  nodes) 
show  much  better  agreement  with  the  exact  solution. 

Finally,  the  solutions  obtained  with  the  finite  element  model  developed  on  this  project 
were  compared  to  those  obtained  with  ANSYS.  Good  agreement  was  also  obtained  in  these 
cases. 

2.3.2  Calculation  of  the  Gradient  of  the  Displacement  Figure-of-Merit 

As  discussed  previously,  given  a  set  of  measured  (or  simulated)  displacement  data,  one 
must  adjust  the  assumed  Young’s  moduli  to  make  the  calculated  displacements  agree  with  the 
measured  values.  This  can  be  formulated  as  an  optimization  problem,  in  which  a  measure  of  the 
disagreement  is  minimized.  A  common  measure  is  the  figure-of-merit  shown  in  Equation  (1). 

Most  efficient  and  robust  techniques  for  minimizing  functions  of  the  form  of  Equation  (1) 
require  information  on  the  change  in  the  predicted  displacements  caused  by  a  small  change  in 
input  parameters.  For  example,  the  Levenberg-Marquardt  technique  requires  extensive 
information  that  is  time-consuming  to  calculate,  i.e.,  the  terms  shown  in  Equation  (2).  The 
BFGS  technique  requires  less  gradient  information,  but  it  is  still  critical  to  develop  this 
information  in  an  efficient  manner.  For  example,  the  crudest  approach  to  calculate  the  gradient 
terms  shown  in  Equation  (3)  would  involve  making  a  small  change  in  a  single  modulus  Ek  and 
then  recalculating  the  entire  displacement  solution.  Repeating  this  process  p  times  to  calculate 
all  such  terms  would  again  require  a  time  on  the  order  of  p  times  n^  or  roughly  n'^.  In  other 
words,  this  brute-force  process  would  require  the  same  length  of  time  as  the  Levenburg- 
Marquardt  method. 

Note  that  the  brute-force  technique  calculates  the  change  in  each  displacement  caused  by 
a  change  in  the  input  Young’s  moduli.  Since  we  then  sum  over  these  displacements  to  calculate 
the  figure-of-merit  J,  we  see  that  the  brute-force  technique  calculates  much  more  information 
than  we  really  need  for  BFGS.  A  better  scheme  is  to  use  an  adjoint  technique  that  efficiently 
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calculates  only  the  needed  information  (Kenton,  1981).  If  we  have  a  set  of  algebraic  equations 
of  the  form: 


0  =  g,(.A:) 


(14) 


together  with  a  single  figure-of-merit  J(x),  one  defines  an  adjoint  x*  associated  with  each 
variable  x  by  constructing  the  following  equations: 


(15) 


After  the  terms  x*  have  been  calculated  using  Equation  (15),  the  terms  of  Equation  (3)  can  be 
obtained  very  simply  by  calculating  the  following  sums: 


_  y  ^Si 

dE,  ' 


(16) 


In  our  case.  Equations  (14)  are  linear,  and  the  matrix  terms  in  Equation  (15)  put  into  the  standard 
form  of  Equation  (9)  are  of  the  form: 

A.  =  km,,  +  0mKM„  - (17) 
ox, 

Since  the  stiffness  and  mass  matrices  KM  and  MM  are  symmetric  (Smith  and  Griffiths,  1998), 
the  matrix  A'  used  to  solve  for  the  adjoints  is  the  same  as  the  matrix  A  used  to  solve  for  the 
displacements. 

This  symmetry  offers  the  potential  for  a  still  larger  savings  in  computation  time,  if 
sufficient  memory  is  available  to  use  direct  sparse  matrix  techniques  rather  than  an  iterative 
scheme  such  as  PCG.  If  we  use  such  a  direct  matrix  solution  technique,  we  will  manipulate  A  in 
some  fashion  (e.g.,  decompose  it  into  upper-  and  lower-diagonal  matrices)  that  will  then  allow 
the  matrix  equation  to  be  readily  solved.  The  time-consuming  part  of  the  solution  process  is  the 
matrix  manipulation  step,  which  takes  a  time  on  the  order  of  the  cube  of  the  rank,  i.e.,  in  this  case 
n^.  The  subsequent  “back-substitution”  process  needed  to  obtain  the  solution  x  for  a  given  b 
takes,  by  contrast,  only  a  time  on  the  order  of  n^.  Thus,  direct  methods  are  very  appealing  when 
used  with  the  BEGS  algorithm:  after  the  matrix  A  has  been  manipulated  to  obtain  a  solution  of 
the  displacements,  the  adjoints  can  be  obtained  with  the  same  matrix  decomposition  in 
essentially  negligible  additional  time.  The  gradient  terms  are  then  obtained  trivially  using 
Equation  (16). 

2.3.3  Determining  Young’s  Modulus  from  Displacements 

The  finite  element  code  was  augmented  to  calculate  the  gradient  terms  after  the 
displacements  were  calculated  using  the  adjoint  method.  The  change  in  figure-of-merit  caused 
by  a  change  in  one  or  more  of  the  Young’s  moduli  calculated  using  brute-force  techniques  was 
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found  to  agree  very  well  with  the  results  of  the  adjoint  method.  The  gradient  information  was 
then  used  to  minimize  the  figure-of-merit  Equation  (1).  We  tried  several  algorithms  to  achieve 
this,  and  the  best  results  were  obtained  using  a  commercially-available  implementation  of  the 
BEGS  algorithm  (IMSL,  1994).  Figure  3  shows  the  relative  rate  at  which  the  figure-of-merit 
converged  in  a  typical  problem  (a  value  of  zero  for  the  figure-of-merit  implies  perfect 
convergence). 

Figure  4  shows  the  Young’s  modulus  obtained  from  computer-generated  MRI 
displacement  data.  Computer-generated  data  were  obtained  by  merely  running  the  finite  element 
code  for  a  specified  Young’s  modulus  distribution,  in  this  case  a  uniform  distribution  with  a 
scaled  value  of  1  except  for  a  square  “tumor”  of  value  10  in  the  upper-right-hand  comer.  The 
minimization  procedure  then  begins  with  a  uniform  (and  incorrect)  Young’s  modulus 
assumption,  and  iterates  until  the  desired  accuracy  is  obtained.  As  can  be  seen  in  Figure  4,  the 
algorithm  does  an  excellent  job  inferring  Young’s  modulus  from  the  computer-generated 
displacement  data.  This  simulation  required  approximately  12  hours  of  computer  time  on  a 
relatively  slow  400  MHz  personal  computer.  Numerous  possibilities  for  increasing  the  speed  of 
this  procedure  are  available  and  will  be  addressed  in  the  next  phase  of  the  project. 

One  detail  of  the  algorithm  needs  to  be  clarified.  Initial  numerical  experiments  resulted 
in  a  much  poorer  reconstmction  of  the  Young’s  modulus  than  is  shown  in  Figure  4.  Inevitably, 
we  found  that  the  method  converged  on  a  solution  that  involved  non-physical  (negative)  Young’s 
moduli  adjacent  to  the  “tumor.”  Further,  the  calculated  Young’s  moduli  inside  the  tumor 
boundary  compensated  for  this  and  were  too  large.  Other  numerical  artifacts  involving  negative 
Young’s  moduli  also  sometimes  appeared  near  the  boundary  of  the  simulated  breast.  To  solve 
these  problems  without  instituting  the  full  complication  of  a  constrained  minimization  approach, 
the  Young’s  moduli  were  mapped  to  a  new  variable  a ,  and  convergence  was  obtained  on  a. 
The  mapping  used: 


(18) 

guarantees  that  the  Young’s  modulus  will  be  larger  than  some  physically-based  minimum  value 
Emin-  This  modification  to  the  algorithm  proved  very  beneficial. 

2.4  Importance  of  Poisson’s  Ratio 

2.4.1  The  Effect  of  Poisson’s  Ratio  on  Compressibility  and  Sound  Speed 

The  example  shown  above  used  prototypical  values  of  the  Young’s  modulus  for  normal 
breast  tissue  and  for  tumors.  However,  the  value  of  Poisson’s  ratio  v  was  0.3,  a  typical  value  for 
structural  (non-biological)  materials.  To  see  what  the  appropriate  values  are  for  biological 
tissues,  it  is  helpful  to  consider  the  speed  of  longitudinal  and  shear  vibrations  in  tissue.  These 
can  be  shown  to  be: 


E{l-v) 
p(l-2v)(l  +  v) 


(19) 
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2/7(1 +  1/) 


(20) 


The  value  for  the  Young’s  modulus  E  of  normal  breast  tissue  is  on  the  order  of  10,000  Pa 
(van  Houten  et  al.  [1999]  use  8000  Pa).  The  density  of  tissue  is  about  that  of  water,  so  p 
-1000  kg/m^.  The  speed  of  longitudinal  vibrations  is  just  the  sound  speed,  which  is  about 
1500  m/sec.  Using  Equation  (19),  we  can  readily  determine  that  Poisson’s  ratio  must  be  very 
close  to  0.5,  so  that  Equation  (19)  simplifies  to: 


c 


2 

L 


E 

3/7(1 -2v) 


(21) 


Equation  (21)  implies  that  v  -  0.499999. 

It  is  interesting  to  note  that  while  the  modulus  E  varies  over  a  reasonably  wide  range,  the 
sound  speed  is  found  to  be  quite  constant,  even  in  tumors  (this  is  why  ultrasound  is  only  of 
limited  use  for  detecting  breast  cancer).  Equation  (21)  thus  implies  that  Poisson’s  ratio  must 
vary  slightly  as  E  varies  to  keep  constant. 

Again  following  Smith  and  Griffiths  (1998),  the  stiffness  matrix  KM  is  developed  from  a 
set  of  constitutive  equations  that  relate  strains  e  to  stresses  o.  In  the  plane  strain  approximation, 
the  associated  matrix  equation  is  written: 


a  =  De 


(22) 


where  D  has  terms  of  the  form 


(23) 


Thus,  as  v/ approaches  0.5,  this  matrix  becomes  ill-conditioned.  Special  techniques  must  be  used 
to  obtain  a  converged  solution,  which  in  any  case  can  be  difficult.  In  the  limit  where  the 
Poisson’s  ratio  becomes  0.5,  the  material  becomes  incompressible,  the  longitudinal  sound  speed 
becomes  infinite,  and  an  entirely  different  formulation  of  the  equations  of  motion  becomes 
necessary.  In  this  formulation,  an  additional  variable  must  be  introduced,  the  hydrostatic 
pressure. 

Most  finite  element  codes  are  incapable  of  handling  the  purely  incompressible  case,  and 
may  have  difficulties  with  the  nearly-incompressible  case  as  well.  Numerical  experiments  with 
the  Creare-built  code  showed  poor  convergence  properties  as  Poisson’s  ratio  was  increased  to 
0.499  and  above.  Presumably  for  the  same  reason,  some  of  the  MRE  work  reported  in  the 
literature  utilized  a  value  of  0.49  (van  Houten,  1999).  At  first  glance,  using  0.49  instead  of 
0.49999  . . .  appears  fairly  reasonable,  until  one  calculates  the  sound  speeds.  At  a  value  of  0.49, 
for  the  Young’s  modulus  and  density  cited  earlier.  Equations  (19)  and  (20)  yield: 

c^«13m/sec  (24) 
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Note  that  the  longitudinal  sound  speed  is  far  less  than  the  actual  value  of  ~1500  m/sec.  Note  also 
from  Equation  (20)  that  the  shear  speed  is  well-behaved  as  v  approaches  0.5. 

At  an  assumed  excitation  frequency  of  200  Hz,  in  the  range  typical  of  MRE  experiments, 
the  corresponding  wavelengths  are  given  by  the  quotient  of  the  sound  speed  and  the  frequency, 
yielding: 


~  6.5cm 

(26) 

~  0.9cm 

(27) 

The  shear  wavelength  is  comparable  to  the  size  of  the  tumors.  As  mentioned  earlier,  achieving 
small  wavelengths  at  the  low  frequencies  resolvable  by  MRI  is  the  advantage  of  using  shear 
waves  in  the  first  place.  However,  the  longitudinal  wavelength  has  been  made  much  shorter  than 
the  actual  value  at  200  Hz  of  about  7  meters.  In  effect,  the  use  of  a  Poisson’s  ratio  of  0.49  has 
made  the  tissue  much  more  compressible  than  it  actually  is. 

2.4.2  Investigative  Calculations  with  ANSYS 

To  investigate  the  implications  of  using  an  artificially  reduced  value  of  Poisson’s  ratio, 
we  performed  a  series  of  calculations  with  ANSYS.  This  code  was  used  rather  than  the  Creare 
code  since  the  former  more  capably  handles  nearly-incompressible  substances.  The  simulated 
“breast”  in  these  calculations  is  a  square  of  dimension  100  mm  x  100  mm,  excited  along  the 
upper  edge  with  a  forced  displacement  whose  amplitude  varies  sinusoidally  from  zero  at  the  left 
and  right  comers  to  a  maximum  value  of  1  mm  in  the  middle.  In  all  cases,  the  excitation 
frequency  is  200  Hz,  the  damping  constant is  .00016  (Gao,  1995),  and  the  density  p  is 
1000  kg/m^.  The  nominal  modulus  E  is  10000  Pa. 

These  various  cases  are  discussed  below. 

Casel:  Nominal  Case 

The  nominal  case  involves  a  homogeneous  material,  i.e.,  a  stiff  inclusion  representing  a 
tumor  has  not  been  included.  The  Poisson’s  ratio  has  been  set  to  a  value  very  close  to  0.5  that 
will  yield  the  proper  longitudinal  sound  speed.  The  “breast”  is  excited  in  the  up-  and  down- 
direction  along  the  upper  edge. 

A  “snapshot”  of  the  calculated  displacement  amplitudes  at  a  particular  time  is  shown  in 
Figure  5.  The  wavelength  of  the  longitudinal  waves  are  much  larger  than  the  10  cm  maximum 
dimension  of  the  breast.  Also,  there  are  no  inclusions  that  could  possibly  cause  mode  conversion 
from  longitudinal  to  shear  waves.  Thus,  no  “waves”  are  seen,  just  positive  amplitude 
displacements  in  the  y-direction. 
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Case  2:  Shear  Excitation 

This  is  the  same  as  Case  1,  except  that  the  breast  is  excited  in  the  x-  (shear)  direction.  In 
this  case,  Figure  6  reveals  shear  waves  with  the  expected  wavelength  of  about  0.9  cm,  that  are 
damped  fairly  quickly  with  depth.  Note  that  the  black  box  shown  in  the  figure  has  no 
significance  in  this  case  (the  box  denotes  the  location  of  the  inclusion  when  one  is  used  in 
Case  4). 

Case  3;  Artificially  Reduced  Poisson’s  Ratio  with  Longitudinal  Excitation 

This  is  the  same  as  Case  1,  i.e.,  there  is  still  no  stiff  inclusion  and  the  breast  is  excited  in 
the  vertical  direction,  but  in  this  case  a  reduced  Poisson’s  ratio  of  0.49  is  used  as  in  some 
previously  published  research.  Figure?  exhibits  longitudinal  waves  with  a  wavelength  of 
roughly  6.5  cm  as  expected  from  Equation  (26).  The  comparison  of  Figure  7  with  Figure  5 
clearly  illustrates  the  deleterious  effects  of  using  too  small  a  value  for  Poisson’s  ratio. 

Case  4:  Correct  Poisson’s  Ratio,  Stiff  Inclusion,  Shear-Excitation 

This  case  is  similar  to  Case  2  in  that  it  employs  the  correct  Poisson’s  ratio.  Also,  in  this 
case  a  stiff  inclusion  is  provided  at  the  location  of  the  dark  box.  The  stiffness  (increased 
Young’s  modulus)  of  the  inclusion  leads  to  an  increase  in  sound  speed.  The  corresponding 
increase  in  shear  wavelength  can  be  seen  in  Figure  8.  The  evident  distortions  in  the  wave 
patterns  provide  the  raw  data  necessary  for  an  inverse  method  to  estimate  the  Young’s  modulus. 

2.4.3  Conclusions  of  Poisson’s  Ratio  Investigations 

The  assumed  Poisson’s  ratio  exhibits  a  profound  effect  on  the  displacement  patterns. 
While  the  use  of  a  reduced  value  (e.g.,  0.49)  greatly  assists  in  numerical  convergence,  this  cannot 
be  recommended  since  it  gives  rise  to  relatively  low  longitudinal  wave  speeds,  correspondingly 
short  wavelengths,  and  substantial  distortion  of  the  calculated  displacements.  An  intriguing 
question  is  whether  this  effect  may  partly  explain  the  artifacts  reported  in  the  literature  when  in 
vivo  data  and  that  from  phantoms  are  used.  To  date,  these  problems  have  been  attributed  to  the 
use  of  plane  strain  approximation  or  linear  materials  models.  If  more  detailed,  3-D  viscoelastic 
models  are  needed,  this  will  further  exacerbate  the  problems  of  long-running  inverse  calculations 
for  any  future  clinical  applications,  so  it  would  be  useful  to  resolve  this  question. 

Since  the  longitudinal  wavelengths  at  the  frequencies  of  interest  are  much  larger  than  the 
breast,  it  appears  entirely  justified  to  use  a  completely  incompressible  formulation  if  this  will 
improve  numerical  performance.  This  should  eliminate  the  numerical  problems  associated  with 
modeling  nearly-incompressible  behavior.  Such  an  approach  was  previously  recommended  by 
Skovoroda  et  al.  (1995). 

2.5  Plans  for  the  Second  Year  of  the  Project 

During  the  second  year  of  the  project,  we  plan  to  implement  an  incompressible 
formulation  of  the  tissue  equations  of  motion  in  the  Creare  code.  We  anticipate  that  this  will 
eliminate  the  numerical  problems  associated  with  the  use  of  a  Poisson’s  ratio  near  0.5.  We  also 
plan  to  undertake  efforts  to  further  improve  the  already  relatively  fast  execution  speed.  Finally, 
we  will  also  attempt  to  obtain  existing  measured  MRI  displacement  data  for  a  phantom  from  a 
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cooperating  institution.  Using  these  data,  we  will  attempt  to  reconstruct  the  Young’s  modulus 
distribution  and  compare  it  to  the  known  values. 

3.0  KEY  RESEARCH  ACCOMPLISHMENTS 

•  An  adjoint  method  developed  on  this  project  provides  an  extremely  efficient  technique 
for  calculating  the  gradient  of  the  goodness-of-fit  measure  that  is  the  basis  of  the  inverse 
algorithm. 

•  Combining  the  adjoint  method  with  the  BFGS  algorithm  provides  an  efficient  means  for 
calculating  tissue  stiffness  from  MRI  displacement  data. 

•  Much  of  the  experience  with  Magnetic  Resonance  Elastography  previously  reported  in 
the  literature  has  utilized  an  incorrect  value  of  Poisson’s  ratio. 

•  Use  of  an  incorrect  Poisson’s  ratio  leads  to  a  drastic  underprediction  of  the  longitudinal 
acoustic  wavelength,  which  in  turn  can  be  expected  to  cause  artifacts  in  the  reconstructed 
hardness  values. 

4.0  REPORTABLE  OUTCOMES 

There  were  no  reportable  outcomes  in  the  first  year  of  the  project. 

5.0  CONCLUSIONS 

Magnetic  Resonance  Elastography  is  potentially  a  very  effective  tool  for  the  early 
diagnosis  of  breast  cancer.  The  lack  of  an  efficient  and  robust  means  for  solving  the  inverse 
problem  is  the  most  challenging  technical  problem  preventing  widespread  clinical  use  of  MRE. 
During  the  first  year  of  the  project,  we  have  developed  and  demonstrated  a  new,  highly  efficient 
algorithm  for  solving  the  inverse  problem.  The  results  obtained  so  far  with  this  algorithm  are 
very  promising.  We  have  also  determined  that  poor  MRE  results  reported  in  the  literature  may 
have  been  partly  caused  by  using  incorrect  values  for  the  Poisson’s  ratio  of  breast  tissue.  During 
the  next  year,  we  will  develop  a  mathematical! y-robust  algorithm  that  will  allow  more  accurate 
modeling  of  tissue  motion  during  MRE. 
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7.0  APPENDICES 


Distance  from  Surface,  m 


Figure  1 .  Comparison  of  Finite  Element  Code  Predictions  of  the  Logarithm  of  the  Displacement 
Amplitude  to  the  Logarithm  of  the  Displacement  Amplitude  for  the  Exact  Solution.  These  data  show  that 
our  finite  element  code  closely  matches  the  exact  solution  for  the  analyzed  configuration. 
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Figure  2.  Instantaneous  Displacements  Versus  Depth  for  Various  Meshing  Schemes.  These  data  show 
the  effect  of  model  order  on  the  match  to  the  exact  solution.  For  this  model,  3600  quadrilateral  elements 
or  20,000  triangular  elements  yield  good  model  accuracy. 
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Reconstructed  Young's  Modulus 
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Figure  4.  Young’s  Modulus  Calculated  from  Computer-Generated  MRI  Displacement  Data  Using  the 
Creare  Code  and  Algorithm.  The  stiffness  (relative  modulus  10)  of  the  simulated  tumor  in  the  upper-right- 
hand  corner  can  be  easily  discerned. 


17 


TM-2136 


€reare 


ANSYS  5.7 
MAY  22  2001 
08:56:26 


NODAL  SOLUTION 

STEP=1 

SUB  =1 

FREQ=200 

UY  (AVG) 

RSYS=0 

PowerGraphics 

EFACET=1 

AVRES=Mat 


DMX 

=.657E-03 

SMN 

=-.854E-04 

SMX 

=.500E-03 

-.854E-04 

PV 

-.510E-04 

-.165E-04 

IMI 

.179E-04 

I^H 

.523E-04 

ir^ 

.867E-04 

.121E-03 

.156E-03 

IFF] 

.190E-03 

.224E-03 

.259E-03 

.293E-03 

.328E-03 

n 

.362E-03 

□ 

.397E-03 

n 

.431E-03 

f  1 

.465E-03 

■1 

.500E-03 

Figures.  Calculated  Longitudinal  Displacements  for  the  Nominal  Case  (longitudinal  excitation)  for 
Correct  Material  Properties  Using  ANSYS. 
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Figure  6.  Calculated  Shear  Displacements  for  Shear  Excitation  with  Correct  Material  Properties  Using 
ANSYS. 
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Figure?.  Calculated  Longitudinal  Displacement  for  Longitudinal  Excitation  with  Artificially  Lowered 
Poisson’s  Ratio  Using  ANSYS. 
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Figure  8.  Calculated  Shear  Displacements  for  Shear  Excitation  with  Correct  Material  Properties  Using 
ANSYS  with  Stiff  Inclusion. 
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